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ABSTRACT 

A  forward  model-based  algorithm  for  reconstructing  images  from  Fourier  domain  measurements  is  described.  The 
model  accurately  treats  broadening  of  the  Fourier  domain  sampling  function  resulting  from  finite  aperture  size  and 
spectral  bandwidth.  The  forward  modeling  approach  also  allows  weighting  of  measurements  according  to  their 
individual  signal-to-noise  ratios  (SNRs),  and  injection  of  a  priori  expectations  about  the  object  such  as  its  support 
and  power  spectral  density.  Image  quality  performance  of  a  large  telescope  array  is  analyzed  as  a  function  of  SNR 
and  array  scale.  The  approach  is  designed  to  enable  joint  processing  of  mixtures  of  intensity  and  amplitude 
interferometric  measurements  and  conventional  incoherent  images,  but  the  focus  of  this  paper  is  limited  to  intensity 
interferometry.  Basic  structures  of  the  object  are  seen  to  be  recovered  as  SNR  approaches  10  for  unity  coherence. 


1.  INTRODUCTION 

Intensity  interferometry  was  proposed  in  the  mid-1950s  by  Hanbury  Brown  and  Twiss  in  the  mid  1950s  [1-3].  Over 
the  ensuing  decades,  the  concept  was  developed  into  an  instrument  that  successfully  measured  the  widths  of  a 
number  of  stars,  and  has  also  been  applied  to  other  fields  outside  of  astronomy  [4-8].  Emphasis  later  shifted  away 
from  intensity  interferometry  as  it  became  apparent  that  amplitude  interferometry  had  better  signal-to-noise  ratio 
(SNR)  and  image  formation  capabilities,  especially  in  view  of  improved  optical  technologies  such  as  adaptive  optics 
[9-11].  In  addition  to  SNR,  image  formation  has  been  a  question  for  intensity  interferometry.  Intrinsic  ambiguities 
were  present  in  reconstructions  of  amplitude-only  data  from  conventional  approaches  such  as  phase  retrieval  [12- 
13].  However,  a  number  of  techniques  were  developed  to  retrieve  an  image  that  worked  under  high-SNR  conditions 
[14-19]. 

In  order  to  address  performance  under  low  SNR  conditions,  new  approaches  have  been  developed  including  forward 
model-based  techniques  [20-22].  A  forward  model  is  a  set  of  mathematical  relations  that  relate  the  pristine  object, 
sensor  design,  and  environmental  conditions  to  the  measurements.  A  forward  model  in  image  reconstruction 
involves  estimating  the  object  that  when  propagated  through  the  forward  model  best  fits  the  data  in  terms  of  some 
cost  function,  subject  to  various  constraints  and  a  priori  expectations,  for  example  non-negativity,  support,  entropy, 
total  variation,  image  sharpness,  and  prior  knowledge.  In  the  most  rigorous  case  the  quality  of  fit  to  the  data  is 
quantified  by  its  likelihood  according  to  a  physics-based  model  of  the  measurement  noise  statistics  characterized  by 
a  joint  probability  density  function  (PDF).  When  this  detailed  knowledge  is  unavailable  or  impractical,  a  weighted 
least  squares  cost  function  is  often  used.  Typically  the  minima  of  the  PDF  cannot  be  solved  for  in  closed  form,  and 
is  therefore  found  by  iterative  search  techniques. 

Forward  models  have  significant  benefits  over  other  reconstruction  techniques,  including  known  optimality 
properties,  ability  to  tolerate  incomplete  and  noisy  data,  and  computation  of  Cramer-Rao  bounds.  They  provide  a 
straightforward  approach  to  factor  the  complexities  of  real  sensors  and  environments  into  the  estimation  process. 

The  primary  disadvantage  of  forward  models  is  their  computational  complexity;  however,  modern  computers  and 
algorithms  have  largely  mitigated  this  disadvantage.  Another  challenge  can  be  occurrence  of  multiple  minima  in  the 
cost  function  which  can  make  finding  the  global  minima  difficult  or  impractical,  in  turn  resulting  in  sub-optimal 
solutions  that  can  depend  on  the  initial  estimates  of  the  free  parameters.  Forward  models  naturally  address  several 
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challenges  in  intensity  interferometry:  complex  non-linear  relation  of  measurements  of  the  object’s  squared  Fourier 
magnitude  to  a  spatial  domain  representation  of  the  object,  sparse  and  irregular  spacing  of  the  measurements  in  the 
Fourier  plane,  and  low  SNR.  Forward  models  also  offer  the  potential  for  computation  of  Cramer-Rao  bounds 
(CRBs)  [23].  CRBs  are  of  special  interest  for  intensity  interferometry  and  other  forms  of  image  formation  because 
they  indicate  fundamental  statistical  limits  to  image  reconstruction  performance  and  can  be  adapted  to  include  the 
effects  of  constraints  and  biases  [23-32]. 

This  paper  develops  a  forward  model  approach  for  image  reconstruction  from  interferometric  data.  The  resulting 
algorithm  is  referred  to  as  FIIRE  (Forward  model  Interferometry  Image  Reconstruction  and  Estimation).  Section  2 
describes  the  forward  model  which  includes  treatment  of  effects  of  finite  aperture  and  spectral  bandwidth  (BW)  on 
broadening  the  Fourier  domain  measurement  response  function  beyond  the  idealized  delta  response.  Section  3 
describes  the  maximum  a  posteriori  (MAP)  cost  function  and  the  iterative  search  approach.  Section  4  describes  a 
conceptual  intensity  interferometry  imaging  system  that  represents  the  extreme  of  current  feasibility  and  Section  5 
characterizes  the  SNR  it  might  achieve  on  a  geosynchronous  satellite.  Section  6  investigates  image  reconstruction 
quality  as  a  function  of  measurement  SNR.  Section  7  summarizes  the  conclusions. 

2.  FORWARD  MODEL 


Intensity  interferometry  operates  by  capturing  light  at  two  or  more  telescopes  and  combining  the  signals.  The  vector 
separation  between  a  pair  of  telescopes,  projected  into  the  plane  orthogonal  to  a  line  between  the  array  and  the 
object,  is  termed  the  baseline.  Light  collected  at  baseline  vector  B  and  wavelength  X  provides  information  about  the 
Fourier  transform  at  spatial  frequency  u  =  B/A.  For  conciseness,  in  the  remainder  of  this  document  the  Fourier 
domain  will  be  referred  to  as  the  UV  plane,  and  a  measurement  of  the  value  of  the  object’s  Fourier  transform  will  be 
referred  to  as  a  UV  sample  or  visibility.  In  amplitude  interferometry  the  light  from  the  telescopes  is  optically 
interfered  allowing  measurement  of  both  the  UV  amplitude  and  phase,  also  referred  to  as  the  complex  visibility  [33]. 
In  phase  closure  techniques,  appropriate  products  of  complex  UV  samples  from  three  or  more  baselines  are  formed 
which  cancel  out  the  effects  of  phase  fluctuations  from  atmospheric  turbulence  and  the  collection  system  [34,35].  In 
intensity  interferometry  the  fluctuations  of  the  intensity  of  the  light  are  measured  at  very  high  bandwidth.  The 
normalized  correlation  of  the  fluctuations  of  the  signals  from  two  telescopes  provides  a  direct  measurement  of  the 
magnitude  squared  of  complex  visibility,  hereafter  referred  to  as  visibility-squared  denoted  \V\2  [2,33]: 


|U|2  (y  _  *k,k'\  _  <AIk(U)AIk,(U)) 
1  1  V  \  )  <ik(U)><ik,(U)> 


(1) 


where  t  is  time,  X  is  wavelength,  telescope  positions  are  indexed  by  k,  I  is  intensity,  A  indicates  fluctuations  about 
the  mean  and  Bk  k,  is  the  baseline  formed  by  telescope  positions  k  and  k\ 

To  make  the  forward  model  as  general  as  possible  the  core  code  does  not  assume  any  particular  geometry  of  the 
telescope  array  configuration.  Its  inputs  are  a  list  of  measurements,  each  specified  by  UV  coordinate  and  noise 
standard  deviation.  The  FIIRE  code  can  accept  any  mixture  of  complex  visibility  and  |  V\2  measurements,  the  set  of 
UV  coordinates  for  each  need  not  match  nor  have  any  regular  spacing,  they  are  completely  arbitrary  for  each  sensing 
modality.  FIIRE  has  been  designed  with  placeholders  to  expand  its  capability  to  utilize  amplitude  only  and  phase 
closure  measurements.  However,  the  focus  of  this  paper  is  limited  to  intensity  interferometry. 

Under  the  conditions  of  negligible  detector  dark  current,  the  measurement  noise  for  a  single  intensity  interferometry 
\V\2  measurement  is  [4], 

I  Pb( ^))  ,7, 

M2  Jf  A/72  p2T(X) 

where  T  is  the  collection  time  over  which  the  signals  are  correlated,  A/  is  the  detector  BW,  pj  and  pB  are  the 
signal  levels  of  the  light  from  the  target  and  any  background  sources  in  photo-electrons  per  second  per  spectral 
bandwidth  (Hz),  and  x  is  the  excess  noise  factor  of  the  detector.  The  photon  rate  spectral  densities  are  proportional 
to  target  and  background  brightness,  telescope  aperture  area,  transmission  of  the  optical  system,  and  detector 
quantum  efficiency.  This  expression  is  accurate  under  conditions  of  low-degeneracy  [33],  which  is  satisfied  by 
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direct  or  reflected  thermal  sources  in  the  optical  regime  except  at  extraordinarily  high  temperatures,  and  detector 
response  BW  much  smaller  than  spectral  BW  of  the  light.  Note  that  the  noise  is  dependent  on  the  spectral  density  of 
the  light  only,  so  once  the  spectral  BW  significantly  exceeds  A/ there  is  no  gain  in  SNR  by  detecting  a  larger  spectral 
BW.  If  multiple  measurements  are  collected  such  that  their  wavelengths  and  baselines  map  to  UV  coordinates  with 
separations  small  compared  to  the  correlation  length  of  the  complex  spectrum,  it  simplifies  things  at  no  disadvantage 
to  average  them  together  to  form  a  single  input  to  image  reconstruction  processing.  The  spectrum  correlation  length 
is  given  byAwcor  ~  1/^Fobj,  where  IFobj  characterizes  the  width  of  the  object.  The  highest  reduction  in  noise  is 
achieved  by  setting  the  weights  equal  to  the  reciprocal  of  the  noise  variances  of  the  individual  measurements,  e.g. 
wm  =  1  /  cr^.  The  effective  noise  standard  deviation  of  the  combined  measurement  is: 

jEjwmam  Jl \vn^l°rn  1 

=  - =  y  1/a2  =  I  (3) 

Lml/Vm  1/cr^ 

If  the  variances  of  all  the  measurements  are  similar  this  results  in  the  standard  VN  type  noise  reduction. 

Now  the  object  spatial  intensity  distribution  needs  to  be  related  to  the  UV  samples.  In  the  ideal  case  a  \V\2 
measurement  samples  a  discrete  point  of  the  magnitude  squared  of  the  object’s  Fourier  transform.  However,  as 
illustrated  in  Fig.  1(a),  two  telescope  apertures  actually  measure  light  over  a  continuum  of  baselines  because  of  their 
finite  size.  The  measurement  is  an  integral  over  the  object’s  \  V\2  spectrum,  with  a  weighting  function  equal  to  the 
autocorrelation  of  the  aperture’s  transmission  function.  For  a  circular  aperture  of  diameter  D ,  the  width  of  the  UV 
footprint  of  the  weighting  function  is  2 D/X.  The  spectral  BW  of  the  signal  causes  a  similar  broadening,  described 
by: 


|AU|  =  Gf - r“)  1*1  =  (Y max  ~  Ymin)\B\  =  AY|S|  (4) 

V/Lmin  Amax/ 

where  the  weighting  function  is  now  a  line  segment  of  length  |Au|  extending  in  a  radial  orientation  from  the  UV 
plane  origin,  and  Ay  is  the  spectral  BW  in  m"1.  For  small  spectral  BWs  and  moderate  apertures  sizes  the  composite 
of  the  two  effects  is  well  approximated  by  their  convolution.  The  resulting  UV  domain  measurement  response 
function,  hereafter  referred  to  as  the  UVPSF,  is  illustrated  in  Fig.  1(b): 

Dj  =  f  du  \0(u)\2  UVPSFj(u)  (5) 

UVPSF Au)  =  f*Jmax  dX  [A  ®  A }  (f)  (6) 

J  A'jmin  VA/ 


where  j  denotes  a  particular  measurement,  A(x )  is  the  transmission  function  of  the  telescope  aperture,  ®  denotes 
correlation,  and  A j>min/max  are  the  spectral  range  of  the  measurement.  An  object’s  UV  spectrum  is  correlated  over 
distances  of  |Au|~l/n  where  H  is  its  angular  subtense.  It  is  desirable  for  the  UVPSF  to  be  much  smaller  than  this 
width. 


a) 


b) 
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Fig.  1.  (a)  Illustration  of  UVPSF  broadening  due  to  finite  aperture  size  of  a  telescope  pair  sensing  light  at  a 
continuum  of  baselines,  (b)  Illustration  of  an  aperture  and  spectrally  broadened  UVPSF  overlaid  on  a  logarithmic 

scale  display  of  an  object’s  UV  spectrum. 

The  FIIRE  forward  model  approximates  the  integral  in  Eq.  (5)  as  follows.  The  object  is  represented  as  a  grid  of 
spatial  domain  values.  It  is  transformed  into  the  UV  domain  by  the  discrete  Fourier  transform  (DFT).  UV  domain 
quantities  are  represented  on  the  rectilinear  grid  of  samples  corresponding  to  object’s  DFT  which  is  hereafter 
referred  to  as  the  UV  grid.  The  values  of  the  UV  PSF  for  each  measurement  j  are  calculated  on  the  UV  grid  using  a 
discrete  approximation  to  Eq.  (6).  Each  UV  PSFy  is  lexicographically  ordered  into  a  row  vector,  and  these  are 
stacked  into  a  single  matrix,  R.  The  size  of  this  matrix  can  be  very  large,  however,  typically  only  a  tiny  fraction  of 
the  elements  of  each  row  are  non-zero.  Using  sparse  matrix  representation  techniques  (a  basic  capability  of 
MatLab)  solves  the  data  storage  issue.  The  forward  model  for  the  noiseless  measurement  is  now  simply  expressed 
as 


g  =  RF\o\2  (7) 

where  6  is  a  lexicographical  ordering  of  the  elements  of  the  spatial  domain  representation  of  the  object  intensity 
distribution,  F  denotes  the  DFT  operation,  and  g  is  the  vector  of  noiseless  measurements.  In  code  it  is  more 
efficient  to  calculate  F  \o\2  by  a  2D  FFT  then  re-order  the  result  into  a  vector.  With  R  pre-calculated,  the  forward 
model  operation  is  extremely  efficient.  The  forward  model  for  amplitude  interferometry  is  nearly  identical  but  with 
|  o  | 2  replaced  by  o  in  Eqs.  (5)  and  (7).  The  case  of  UVPSFs  that  are  smaller  than  UV  grid  spacing  is  currently 
treated  by  setting  the  values  of  R  to  be  equivalent  to  bilinear  interpolation.  The  code  contains  a  placeholder  to 
switch  to  sine  interpolation,  which  by  the  Fourier  sampling  theory  is  rigorously  consistent  with  the  physics. 
Truncating  the  sine  function  to  moderate  size  remains  accurate  yet  still  allows  R  to  be  highly  sparse. 


3.  COST  FUNCTION  AND  SEARCH  ALGORITHM 


The  cost  function  is  given  in  Eq.  (8).  It  consists  of  a  data  agreement  metric,  shown  as  the  first  term,  followed  by 
five  regularizing  terms. 


E  =  YJj^\\Mj\P -\gj{6{x))\P\  -  KEntZx6(x)log(d(x))  -  Knk  £*|6(*)l 


1 2 


1/2 

+  KTV'Lx  (|Vxo(x)I2+  |Vyo(x)|2)  +  KPSD_1D  £n 


PSD( \u\n) 


-D 


(8) 


+  KX 


PSD-2D 


Zu 


\8m\ 


PSD(\u\n) 


D 


In  Eq.  (8),  o(x)  is  the  estimate  of  the  object  and  0(x)  its  DFT,  Mj  denotes  a  noisy  measured  value  and  gj{d(x))  its 
prediction  based  on  the  object  estimate  and  forward  model,  derivatives  of  the  object  in  the  x  and  y  directions  are 
indicated  by  Vx  and  Vy.  The  agreement  metric  compares  the  discrepancy  between  the  measurements  and  their 
predicted  values.  The  discrepancy  is  quantified  as  an  L -a  norm,  where  a  is  between  1  and  2.  There  is  also  a 
shaping  exponent,  /3 ,  which  can  increase  sensitivity  to  small  values  of  the  measured  data.  The  measurement  terms 
are  individually  weighted  according  to  their  respective  variances  of.  As  described  in  more  detail  below,  the 
constant  D  in  the  final  two  terms  controls  whether  they  act  purely  as  smoothness  constraints  by  penalizing  non-zero 
Fourier  values,  or  act  to  encourage  reconstructions  to  be  consistent  with  a  priori  expectations  about  the  object’s 
power  spectral  density  (PSD). 

The  regularizing  terms  in  Eq.  (8)  each  serve  different  functions.  The  entropy  term,  with  leading  factor  KEnh  acts  to 
smooth  the  object  and  encourage  the  energy  of  the  object  to  be  contained  in  a  smaller  number  of  pixels  or 
equivalently  discouraging  small  amounts  of  energy  in  many  pixels.  The  Tikhonov  or  image  sharpness  regularization 
term,  with  leading  factor  KTik,  tends  to  concentrate  energy  in  a  limited  number  of  pixels.  The  next  term  in  Eq.  (8)  is 
the  total  variation  term,  with  leading  factor  KTV.  This  well  known  regularizing  approach  serves  to  suppress  noise 
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while  preserving  sharp  edges.  There  is  not  currently  any  provision  for  optimizing  the  relative  weights  of  each  term, 
this  is  left  to  the  user.  The  term  with  leading  factor  KPSd-id ,  imposes  prior  information  in  terms  of  the  expected  PSD 
as  a  function  of  radial  distance  of  a  UV  coordinate  from  the  origin,  e.g.  \u\.  This  is  accomplished  by  dividing  the 
UV  domain  into  a  series  of  concentric  annuli,  with  Sn  specifying  the  set  of  Nn  UV  grid  elements  in  annulus  n ,  and 
\u\n  characterizing  the  annulus’  radius.  The  average  of  the  object  estimate’s  PSD  values  within  an  annulus  is 
calculated  and  divided  by  the  expected  PSD.  When  D  takes  the  value  1,  this  term  encourages  the  estimated  object  to 
be  consistent  with  expectations  about  its  power  spectrum.  When  D  equals  0,  this  term  encourages  object 
smoothness  by  incurring  a  penalty  for  larger  PSD  values,  with  the  penalty  at  a  spatial  frequency  weighted  by  the 
inverse  of  the  a  priori  PSD.  Since  the  PSD  in  the  denominator  is  the  average  a  priori  PSD  as  function  of  \u\n  it 
does  not  express  any  assumptions  about  the  object’s  rotational  orientation. 

This  cost  function  can  be  extended  to  incorporate  data  from  other  sensing  modalities  such  as  amplitude 
interferometry,  phase  closure,  and  conventional  imaging  by  adding  additional  data  agreement  cost  terms  for  each 
modality  analogous  to  the  1st  term  in  Eq.  (8).  This  enables  joint  processing  of  all  data  to  find  an  object  solution  that 
is  in  balanced  agreement  with  all  data  sources  and  the  regularizing  priors,  while  also  subject  to  non-negative  and  an 
optional  support  constraint. 

A  solution  for  the  grid  of  object  intensity  values  that  minimizes  the  cost  function  is  sought  using  the  L-BFGS-B 
iterative  limited  memory  quasi-Newton  search  algorithm  [36].  As  it  searches  through  parameter  space  L-BFGS-B 
accelerates  the  search  by  developing  an  approximate  model  of  the  local  2nd  order  partial  derivatives  without 
requiring  evaluation  or  storage  of  the  full  Hessian  matrix.  It  also  has  provisions  for  placing  upper  and  lower  bounds 
on  each  parameter  providing  a  direct  approach  to  enforcing  the  non-negative  property  of  the  object.  The  FIIRE  code 
allows  this  option  or  alternatively  a  re-parameterization  of  the  object  as  the  square  of  an  underlying  function  which 
also  has  the  effect  of  ensuring  non-negativity.  Support  constraints  are  easily  employed  by  limiting  the  set  of  object 
grid  elements  which  are  optimized  and  leaving  the  rest  at  zero.  This  approach  requires  calculations  of  the  gradients 
of  the  cost  function  with  respect  to  the  object;  closed  form  expressions  for  each  term  in  Eq.  (8)  are  easily  calculated 
but  not  shown  here.  Use  of  Eq.  (7)  allows  the  efficient  simultaneous  computation  of  all  elements  of  the  gradient  of 
the  1st  term  of  Eq.  (8),  the  data  agreement  term,  using  FFTs  and  matrix  multiplies.  One  other  regularization 
technique  supported  by  FIIRE  that  should  be  mentioned  is  to  re-parameterization  the  object  as  the  convolution  of  a 
low-pass  kernel  with  an  underlying  function.  This  has  been  referred  to  as  the  method  of  sieves  [38]  and  effectively 
places  a  support  constraint  on  the  object  estimate’s  UV  spectrum  preventing  amplification  of  high-frequency  noise. 
If  the  transfer  function  corresponding  to  the  Fourier  transform  of  the  kernel  decreases  with  increasing  \u\  then  it  will 
act  to  attenuate  the  convergence  of  higher  frequency  components  which  can  help  to  avoid  local  minima.  It  should 
be  noted  that  the  algorithmic  approach  used  by  FIIRE  is  quite  similar  to  that  of  MIRA  [22,  34],  another  image 
reconstruction  code  widely  recognized  in  the  imaging  interferometry  community  [38]. 


4.  INTENSITY  INTERFEROMETRY  SYSTEM 

To  demonstrate  the  FIIRE  code  and  investigate  influence  of  its  various  features  on  reconstructed  image  quality,  an 
intensity  interferometry  collection  system  was  defined  based  on  the  Cherenkov  Telescope  Array  (CTA)  which  is 
comprised  of  97  telescopes  arranged  on  a  1  km  wide  square.  This  system  is  described  in  detail  in  reference  [39]. 
Note  the  pattern  of  telescope  positions  is  chosen  to  give  a  higher  diversity  of  baselines  than  would  a  pattern  with 
equal  spacing.  Nevertheless,  many  of  the  telescope  pairs  have  equivalent  baselines  resulting  in  considerable 
redundancy  in  the  UV  measurements  especially  at  smaller  separations.  Since  an  emphasis  of  this  study  was  to 
investigate  ground-based  imaging  of  objects  in  geosynchronous  Earth  orbit  (GEO),  the  square  pattern  of  telescope 
positions  was  scaled  down  to  40  m  across  which  is  better  matched  to  the  angular  subtense  and  desired  angular 
resolution.  Fig.  2  shows  the  geometry  of  this  array  and  the  resulting  baselines  assuming  a  zenith  line-of-sight. 
Color-coding  indicates  the  level  of  redundancy  at  each  baseline.  The  main  details  of  the  system  parameters  are 
summarized  in  the  following  list: 

•  97  telescopes  arranged  within  a  40x40  m  square  providing  653  unique  baselines 

•  0.75  m  diameter  apertures  with  0.375  m  diameter  central  obscurations  for  .3313  m2  collection  area 

•  0.196  net  system  transmission  (atmosphere  *  optics  *  quantum  efficiency) 

•  Nyquist  sampling  for  maximum  horizontal  baseline  is  5  nrad  (7.3E-4  arcsec),  corresponding  to  17.9  cm  at 
GEO 
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•  Light  is  detected  in  84  spectral  channels  spanning  400-900  nm 

•  The  channels  are  evenly  spaced  in  wavenumbers  at  0.0167  pm"1  and  each  is  0.0083  pm'1  wide.  This 
comprises  50%  utilization  of  the  full  400-900  nm  spectrum. 

•  Each  channel  is  divided  into  25  contiguous  sub-channels  of  100  GHz  width  and  spacing 

•  Light  in  each  sub-channel  is  measured  with  a  500  MHz  BW  detector  and  1  GHz  sample  rate 

Figure  3  diagrams  the  division  of  the  broad  spectrum  into  channels  and  sub-channels.  Physical  methods  of 
achieving  this  were  not  a  focus  of  this  investigation,  but  a  potential  approach  that  could  accurately  separate  the 
channels  at  the  necessary  spectral  resolutions  is  use  of  Echelle  cross-grating  spectrometers  combined  with  micro¬ 
lens  assemblies  to  collect  and  couple  the  light  into  optical  fibers  fed  into  high  BW  detectors.  The  goal  of  dividing 
up  the  spectrum  is  to  gain  best  advantage  of  the  full  visible  band  of  light.  Consider  a  large  satellite  of  40  m 
maximum  width  at  a  range  of  40,000  km.  This  corresponds  to  an  angular  subtense  of  1  pxad.  Nyquist  sampling  of 
its  UV  spectrum  is  Au  =  1  pxad'1.  By  Eq.  (1),  at  the  maximum  baseline  of  40  m,  a  spectral  channel  of  0.048  pm'1 
BW  maps  to  a  line-segment  on  the  UV  plane  of  length  0.27  qrad"1  thus  is  significantly  smaller  than  the  correlation 
length  of  the  object’s  UV  spectrum  and  can  be  treated  as  nearly  a  point  sample  such  that  UVPSF  effects  are  not  a 
concern.  On  the  other  hand,  at  400  nm  the  0.75  m  aperture  diameter  results  in  a  UV  PSF  width  of  3.75  qrad"1  and 
could  cause  problematic  blur  of  the  UV  measurements.  To  get  higher  utility  of  the  light  within  a  channel  it  is 
divided  into  25  contiguous  sub-channels  that  are  separately  measured.  Each  pair  of  matching  sub-channels  collected 
by  a  pair  of  telescopes  is  correlated,  then  each  set  of  \V\2  measurements  within  a  channel  are  averaged  to  boost  the 
SNR  by  V25.  This  partially  compensates  for  the  restriction  that  SNR  is  limited  by  the  BW  of  the  detector  response 
not  the  BW  of  the  detected  light.  Fig.  4  shows  the  resulting  set  of  UV  sample  locations  of  the  described  system 
color-coded  by  their  wavenumber. 
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Fig.  2.  Telescope  array  geometry  and  corresponding  set  of  baselines  for  the  analyzed  intensity  interferometer 

system. 


fFor  clarity  of  display  the  number  of 
channels  and  their  widths  has  been 
adjusted  from  actual  values. 
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Fig.  3.  Division  of  broadband  light  into  channels  and  sub-channels. 
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UV  Sample  Locations 


653  unique  baselines 
84  spectral  channels 
54,852  UV  samples  map  to  4905 
distinct  object  UV  correlation  cells 
Compare  to  40,000  pixels  for 
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Fig.  4.  Set  of  UV  samples  collected  by  the  described  system  color-coded  by  wavelength. 

It  is  worth  some  comments  about  the  challenges  for  the  discussed  system.  By  dividing  the  spectrum  into  84 
channels,  each  comprised  of  25  sub-channels,  each  of  the  97  telescopes  would  require  2100  high  BW  detectors,  thus 
203,700  for  the  full  system.  To  ensure  the  spectral  BW  of  matching  sub-channels  from  different  telescopes  overlap 
by  at  least  90%  will  require  exquisite  spectral  accuracy  of  10  GHz  which  corresponds  to  5  picometers  at  the  400  nm 
end  of  the  spectrum.  With  4656  telescope  pairings,  to  avoid  storing  the  203E+3  giga-samples/sec  of  data  that  is 
generated  will  require  approximately  10  million  correlators.  This  design  was  chosen  to  be  representative  of  what  at 
the  time  of  this  paper  is  an  achievable,  but  large,  complex,  and  costly  system. 

5.  Visibility  SNR 

Photo-electron  signal  level  is  now  considered.  The  photon  flux  at  the  top  of  the  atmosphere  for  a  Mv  10  object  is 
illustrated  in  the  left-hand  plot  of  Fig.  5  as  a  function  of  wavelength  [40].  Assuming  a  clear-aperture  area  of  0.3313 
m2  and  net  transmission  efficiency  of  19.6%  including  atmosphere,  optics,  and  quantum  efficiency,  results  in  a 
spectral  density  of  about  7E-10  photo-electrons  per  second  per  Hz.  Photo-electron  rate  spectral  density  is  plotted  vs. 
wavelength  in  the  right-hand  plot  of  Fig.  5. 
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Fig.  5.  Top-of-atmosphere  photon  flux  from  a  Mv  10  object  (left),  and  resulting  photo-electron  spectral  density  at 

the  detector  (right). 


An  individual  measurement  is  defined  as  the  correlation  between  two  sub-channels  which  for  a  1  hour  integration 
time  and  signal  level  of  7E-10  pe-/sec/Hz  results  in  a  |fj2  noise  standard  deviation  of  -1500.  It  simplifies  the 
processing  at  no  loss  to  combine  individual  measurements  that  are  close  together  relative  to  the  UV  spectrum’s 
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correlation  length  by  averaging  them  weighted  by  the  inverse  of  their  variances.  This  was  accomplished  in  this  study 
by  grouping  individual  UV  measurements  into  sets  such  that  all  members  in  a  set  lied  within  0.375  *  A uCOY  of  each 
other.  The  SNR  of  each  resulting  composite  UV  sample  j  is  defined 


SNR\v\2 


m)\2 

aCJ 


(9) 


where  ac j  is  defined  by  Eqs.  (2-3).  The  average  number  of  individual  measurements  per  composite  UV  sample  was 
-2000  resulting  in  an  average  value  of  oc  of  77  which  is  a  factor  of  20  improvement  compared  to  the  noise  of  the 
individual  constituent  samples.  The  noise  level  oc j  and  SNR  will  vary  for  different  composite  UV  samples  Uj  for  a 
number  of  reasons,  (a)  variations  in  the  level  of  redundant  telescope  pairings  with  same  baseline,  (b)  variations  in 
the  number  of  baseline  and  wavelength  combinations  mapping  to  a  particular  Uj,  (c)  variations  of  photon  flux  at 

different  wavelengths,  (d)  variations  in  the  object’s  UV  spectrum  |0(u)|2  which  for  most  objects  will  tend  to 
follow  a  power-law  type  decrease  \u\~a.  Panels  (b-c)  of  Fig.  6  show  a  histogram  of  the  values  of  oc j  and 
characterize  the  variations  with  \uj\.  By  Eq.  (2)  noise  standard  deviation  scales  inversely  to  the  square-root  of  the 
collection  time. 


To  further  characterize  SNR  it  is  useful  to  consider  a  specific  object.  The  pristine  intensity  distribution  shown  in 
Fig.  7(a)  was  generated  from  a  3D  model  and  is  placed  at  a  range  of  35.8E+3  km.  This  object  will  also  be  used  for 
image  reconstruction  analysis  in  Section  7.  Fig.  7  also  shows  (b)  the  object’s  PSD,  (c)  the  set  of  UV  sample 
locations  measured  by  the  system  described  in  Section  2,  and  (c)  the  image  collected  by  a  conventional  full-aperture 
telescope  with  diameter  such  that  its  modulation  transfer  function  (MTF)  cutoff  matches  the  UV  sample  with  the 
largest  radius  indicated  by  the  red  circle  in  (d).  Fig.  6  provides  several  plots  characterizing  the  distributions  of  \V\2, 
noise,  and  SNR,  and  their  trends  with  UV  coordinate  radius  \u\  for  a  1-hour  signal  collection  time.  Per  Eq.  (2)  SNR 
will  scale  in  proportion  to  the  square-root  of  the  collection  time.  It  is  observed  that  noise  level  tends  to  be  higher  for 
larger  \u\  which  results  from  the  array  geometry  providing  fewer  redundant  telescope  pairs  at  longer  baselines. 
Combined  with  a  rapid  decrease  in  |fj2  at  higher  spatial  frequencies,  SNR  drops  rapidly  with  increasing  \u\.  Fig.  6(f) 
indicates  that  for  a  1-hour  collection  time  the  bulk  of  the  UV  samples  have  SNRs  between  10"6  to  10"4. 
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Fig.  6.  Plots  characterize  \  V\2  magnitudes,  noise,  and  SNR  levels  and  their  trends  vs.  radial  distance  from  UV  origin 
based  on  a  1-hour  signal  collection  time  for  the  satellite  object  of  Fig.  7  at  Mv=10  brightness. 
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*  displayed  raised  to  0.5  power 
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Fig.  7.  (a)  Pristine  object;  (b)  object’s  PSD;  (c)  diffraction  limited  image  from  full  aperture  telescope  with  MTF 
cutoff  matching  highest  UV  sample;  (d)  UV  sample  positions. 


6.  IMAGE  RECONSTRUCTION  ANALYSIS 


6.1  Image  Quality  vs.  SNR 

This  section  examines  reconstructed  image  quality  dependence  on  SNR  of  the  |fj2  measurements.  The  test  object 
and  interferometry  system  are  as  described  in  Sections  4  and  5.  Reconstructions  were  performed  using  the  FIIRE 
code.  The  simulation  parameters  are  summarized  in  the  following  list: 

•  The  object  was  estimated  ona256><256  grid  corresponding  to  a  40x40  meter  FOV,  exactly  matching  the 
grid  the  original  pristine  object  used  to  simulate  the  |E|2  measurements  is  represented  on 

•  Iterations  were  performed  until  convergence  to  numerical  tolerance  or  a  maximum  8,000  iterations 

•  The  UVPSD  was  set  to  a  delta  function,  i.e.,  spectral  and  aperture  induced  broadening  were  not  modeled 

•  Initial  estimate  was  a  Gaussian  spatial  distribution  multiplied  by  white  noise 

•  Object  non-negativity  was  enforced  but  a  support  constraint  was  not  used 

•  Earth  rotation  was  ignored,  so  the  orientation  of  the  telescope  array  relative  to  the  object  was  treated  as  fixed 
over  the  observation  time 

•  Images  were  generated  using  each  of  the  different  regularization  terms  activated  one  at  a  time. 

Combinations  were  not  explored.  For  each  term  the  regularization  weight  was  scanned  over  a  large  range  of 
values  to  determine  an  optimum  strength.  The  ID  PSD  term  with  D=  1  (i.e.,  discrepancy  not  penalty)  was 
found  to  produce  the  best  results  (visually).  A  weight  of  IE-3  was  found  to  be  optimum  for  all  cases  except 
infinite  SNR  at  which  point  a  weight  of  IE-5  was  better. 
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Figure  8  shows  the  image  reconstructions  obtained  for  a  range  of  measurement  SNR  levels.  Note  that  lack  of  phase 
information  results  in  an  ambiguity  such  that  a  solution  and  its  180°  rotation  are  equivalent.  For  purposes  of  display 
the  reconstructions  were  rotated  to  match  the  pristine  object  as  necessary.  Two  SNR  metrics  are  annotated  on  each 
image.  SNR|V|=i  is  calculated  from  Eqs.  (2-3)  using  the  average  of  the  25  sub-channels  provided  by  a  single  channel 
from  a  single  telescope  pair  with  |v|2=l,  e.g.  the  object’s  normalized  PSD  value  at  the  origin  of  the  UV  plane,  or 
equivalently  unity  coherence  between  the  intensity  fluctuations.  For  a  1  hour  integration  time,  SNR|V|=i=0.0033. 
The  peak  SNR  for  the  UV  samples  is  characterized  by  SNRDCp  =  88*SNR|V|=i  and  corresponds  to  the  SNR  at  |V|=1 
for  the  composite  UV  sample  with  lowest  noise  standard  deviation.  It  is  observed  that  the  basic  dumbbell  structure 
of  the  object  starts  to  become  discernible  at  SNRDCp  =  1.6.  More  detailed  features  start  to  appear  as  SNRDCp  exceeds 
5.  Near  diffraction-limited  quality  is  achieved  at  SNRDCp  -  29.  Figure  9  is  similar  to  Fig.  8  but  with  the 
reconstructions  performed  without  any  regularization.  It  is  obvious  that  regularization  is  effective  in  suppressing 
meaningless  high-frequency  structure  to  produce  an  image  with  more  pleasing  appearance.  It  is  debatable  whether 
they  are  better  from  the  standpoint  of  the  information  visually  conveyed,  but  the  suppression  of  false  high-frequency 
structure  may  likely  benefit  automated  information  extraction  methods. 


*For  improved  clarity  images  are  displayed  with  dynamic  range  thresholding  and  gamma  adjustment 


Fig.  8.  Image  reconstruction  results  for  a  range  of  SNR  levels. 
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*  For  improved  clarity  images  are  displayed  with  dynamic  range  thresholding  and  gamma  adjustment 


Fig.  9.  Same  as  Fig.  8  but  without  regularization. 

6.2  Optimal  Telescope  Array  Scaling 

Figure  10  compares  reconstructions  obtained  with  the  telescope  location  pattern  scaled  to  1,  V2,  and  %  that  of  Fig.  2. 
Telescope  aperture  area  was  held  fixed.  The  middle  reconstruction  corresponding  to  lA  scale  appears  the  best.  The 
rightmost  is  clearly  worse  and  though  the  left  reconstruction  contains  finer  structures,  they  do  not  correlate  with 
fine-scale  features  in  the  pristine  image  and  are  just  amplified  noise.  This  suggests  there  is  benefit  to  optimizing  the 
array  pattern  to  match  the  characteristic  width  of  the  object,  such  that  the  density  of  samples  and  levels  of  SNR- 
boosting  redundancy  are  concentrated  where  they  are  most  useful. 
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Fig.  10.  Image  reconstructions  vs.  different  scaling  of  the  telescope  array  pattern. 
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6.3  UV  Sample  Weighting 


Three  choices  for  specifying  the  UV  sample  variances,  of,  used  in  the  1st  term  of  the  cost  function  expressed  by  Eq. 
(8)  were  explored.  The  first  was  to  use  the  true  variances  of  the  noise  added  in  simulating  the  visibility 
measurement  data  defined  by  Eqs.  (2-3),  which  matches  the  measurement  physics  and  method  of  combining 
redundant  samples  given  by  Eq.  (3).  The  second  method  was  to  scale  these  by  a  multiplicative  factor.  The  third  was 
to  use  an  identical  value  for  all  UV  samples.  The  third  method  led  to  much  worse  performance  than  the  other 
methods.  Variations  in  image  reconstruction  quality  were  imperceptible  for  scale  factors  between  0.001  -  1000 
beyond  which  quality  began  to  degrade. 


7.  SUMMARY 

A  code  was  developed  for  reconstructing  images  from  long-baseline  interferometry  data  based  on  a  forward  model 
approach  in  which  an  object  is  estimated  such  that  when  propagated  through  a  physics-based  forward  model  of  the 
sensing  system  maximizes  the  match  to  the  set  of  measurements.  The  code  currently  supports  joint  processing  of 
complex  amplitude  and  \  V\2  type  measurements.  Explorations  of  different  regularization  methods  and  measurement 
weighting  schemes  gave  insights  into  the  approaches  that  were  reliably  best.  It  was  observed  that  scaling  the 
telescope  array  location  pattern  to  match  the  resulting  UV  measurement  location  and  SNR  distribution  to  the  PSD  of 
the  object  was  important  to  maximizing  image  quality.  The  performance  of  a  conceptual  system  representing  the 
extreme  of  what  might  be  achieved  with  current  technology  was  considered.  The  system  consisted  of  an  array  of  97 
telescopes,  each  with  a  0.75  m  diameter  aperture  using  2100  narrowband  spectral  channels  spread  across  400-900 
nm  to  make  highly  efficient  use  of  the  photons  within  the  visible  band.  Signals  from  all  4656  pairings  of  telescopes 
were  correlated  to  gain  the  full  Fellgetf  s  advantage  of  the  system.  The  array  pattern  and  wide  spectral  bandwidth 
resulted  in  many  measurements  of  nearly  identical  UV  locations,  which  were  averaged  using  a  scheme  that  resulted 
in  a  typical  number  of  many  hundred  redundant  measurements  per  UV  locations.  Taking  the  SNR  benefit  from  this 
redundancy  into  account,  for  a  Mv  10  object  and  UV  location  with  \V\2  ~1,  it  was  determined  that  an  SNRDCp  of 
approximately  2  was  needed  to  recover  the  coarse  structure  of  the  object.  With  the  described  system  this  would  take 
at  least  1000  hours  of  integration  time.  Finer  details  of  the  object  started  to  appear  as  SNRDCp  exceeded  5,  with 
image  quality  increasing  steadily  up  to  SNRDCr~30  at  which  point  a  near  diffraction-limited  quality  image  was 
obtained.  This  study  suggests  that  in  principle  intensity  interferometry  could  produce  images  of  GEO  objects. 
However  at  typical  GEO  object  brightness  levels,  major  advances  in  measurement  technology  or  algorithmic 
approaches  will  be  needed  to  achieve  images  with  reasonably  small  collection  times.  It  will  be  interesting  to  see  how 
this  picture  changes  with  use  of  other  interferometric  techniques  such  as  amplitude  interferometry. 
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